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amoeboid organism of Physarum plasmodium 
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The amoeboid organism, the plasmodium of Physarum polycephalum, behaves on the basis of 
spatio-temporal pattern formation by local contraction-oscillators. This biological system can be re- 
garded as a reaction-diffusion system which has spatial interaction by active flow of protoplasmic sol 
in the cell. Paying attention to the physiological evidence that the flow is determined by contraction 
pattern in the plasmodium, a reaction-diffusion system having self-determined flow arises. Such a 
coupling of reaction-diffusion-advection is a characteristic of the biological system, and is expected 
to relate with control mechanism of amoeboid behaviours. Hence, we have studied effects of the 
self-determined flow on pattern formation of simple reaction-diffusion systems. By weakly nonlinear 
' analysis near a trivial solution, the envelope dynamics follows the complex Ginzburg-Landau type 

equation just after bifurcation occurs at finite wave number. The flow term affects the nonlinear 
term of the equation through the critical wave number squared. Contrary to this, wave number 
isn't explicitly effective with lack of flow or constant flow. Thus, spatial size of pattern is especially 
important for regulating pattern formation in the plasmodium. On the other hand, the flow term 
is negligible in the vicinity of bifurcation at inflnitely small wave number, and therefore the pattern 
QQ I formation by simple reaction-diffusion will also hold. A physiological role of pattern formation as 

above is discussed. 
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(3 ; I. INTRODUCTION 

vn ' 

, The Plasmodium of Physarum polycephalum is a large amoeboid cell, showing contraction-relaxation cycles every- 
Q.^ ■ where within an organism. These local contraction motivates intracellular transport of endoplasmic sol [Q. Therefore, 
the Plasmodium can crawl when the local flow is appropriately arranged all through the organism. Some types of 
Q • spatio-temporal pattern of the contraction have been observed after stimulation, and discussed in relation to devel- 
C/3 opment of its amoeboid behaviours . Mechanism of the pattern formation is closely related to a mechanism 

of cellular behaviours in the plasmodium. Then we expect that a certain origin of biological fineness for controlling a 
' system is in characteristic mechanism of the biological pattern formation. This paper is concerned with such character 
O |. in a biological system. 

1!^ ' Let us consider theoretical framework of the pattern formation in the plasmodium. Chemical oscillation is a clock 
of the rhythmic contraction. Possible candidates for the chemicals are, for example, Ca^+ and/or ATP Since 
the contraction apparatus is located at outer layer of the plasmodium, called ectoplasm [BIbU1(|, it is proper to focus 
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?H ^ on chemical oscillations and diffusion in the ectoplasm. As above, the plasmodial contraction system can be regarded 
. . . I as diffusively coupled oscillators or a reaction-diffusion system in the ectoplasm. Actually, some experimental results 
can be explained using this framework ^1 12 1. 

But another effect of spatial interaction must be considered. That is the protoplasmic streaming in inner part 
of the Plasmodium, called endoplasm, because the contraction pattern is modified by inhibition of the streaming 
P, p"3p4| ]. Miyake et al. jl^] proposed a model of the information processing system with two levels of subsystems 
corresponding to the endoplasmic oscillators with long-range interaction and the ectoplasmic ones with short-range 
interaction. From more physical viewpoints, some models based on hydrodynamics and chemical kinetics, have been 
presented ||l^,|l^,|l^ . However, we still have an issue that the streaming can be sometimes negligible, as mentioned in 
the former paragraph. To solve the pretence of this discrepancy, it is inevitable to study pattern formation including 
the streaming effect. For that purpose, we firstly suppose a physical model of reaction-diffusion-advection equations. 

The oscillating chemicals are exchanged between endo- and ectoplasm, and flow out/in via the streaming ||l8|,|l9|. 
The streaming is motivated by the gradient of contraction force pO| ]. Hence, the contraction pattern determines 
the streaming. That is to say, the protoplasmic streaming is a factor of the contraction pattern, and also inversely 
motivated by the pattern. Then the reaction-diffusion system having self-determined flow arises. The self-determined 
flow is just a characteristic of the biological system in the plasmodium. We analyze pattern formation of such a system 
by weakly nonlinear analysis, particularly noticing effects of the advection term on simple reaction-diffusion systems. 
The weakly nonlinear analysis is based on singular perturbation expansion, and therefore has an advantage that 
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model equations are free from details. Because, no matter how complex or uncertain the functions are in our model 
equations, expanded lower order terms are merely considered in the analysis. This advantage is valuable especially in 
analysis of biological model, which often includes uncertain assumptions in detail. 

As the results of the analysis, a possible solution of the pretence discrepancy described earlier is found, and 
furthermore, one apparent character of the plasmodial pattern formation is shown. A possible role of the plasmodial 
pattern formation is discussed at a physiological point of view. 



II. BASIC EQUATIONS 

The Plasmodium of Physarum has a cytoplasmic cortex (ectoplasmic gel) filled with endoplasmic sol. The ectoplasm 
makes periodic contraction and relaxation, and it causes intracellular streaming of the endoplasm. A sheet of cytoplasm 
become thick when the endoplasm flowing into it. Metabolic substances regulate contractile cycles, diffuse in the 
cytoplasm, and are transferred by the endoplasmic flow. Standing on these points, we present a dynamics of rhythmic 
pattern formation in the Physarum plasmodium as the following equations for metabolic elements in the ectoplasmic 
gel, Ugci and endoplasmic sol, Ugoi- 



du 



■gol 



Fgel(^u) + V- (i^gclVUgel), (1) 



dt 

— ^ + u ■ Vusoi = Fsoi(/i, u) + - V • {DsoiVhu, 
at h 



.o\) 



Where h denotes the thickness of the endoplasmic sol, and v is the averaged velocity of the endoplasmic flow, u 
is N metabolic species in gel and sol, (Ugci,Usoi) = Fgd and Fgoi represent reaction kinetics among 

metabolic elements and exchanges of them between the ectoplasm and endoplasm. I?goi and Ugoi are diagonal matrices 
of diffusion constants of metabolic elements. We note that Eq. (|l|) becomes the closed system if the dynamics of v, the 
equation of motion of endoplasmic flow, are given. Let us consider the small Reynolds number flow in a contractile 
vein and assume the intracellular pressure is determined by the concentration of the chemicals: v = —q{h)\IP{\i). 
This means the stationary flow approximation that deformation of the ectoplasm is very slow compared with variation 
of the endoplasmic flow. It is hence not considered that the ectoplasmic cortex has viscoelastic features. 

In Eq. p), the sol-gel conversion is ignored, and thus the mass conservation of endoplasmic sol is satisfied. This 
implies the limitation of the model, such as the cell motility, the formation and reconnection of a network of pro- 
toplasmic strands. In the following, we assume that the thickness of the endoplasm is almost constant all over the 
Plasmodium, and that the diffusion constants of the metabolic elements are homogeneous in the plasmodium. The 
intracellular pressure P is expanded around the homogeneous static state u = as 

P(u) — P(us) + y^l"'' - ul)——(us) -f higher order terms. 
^ — ' ou^ 
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Hereafter, we ignore the higher order terms which have no effects on results deduced by the weakly nonlinear analysis. 
Under these assumptions, we rewrite Eq. in the form of the rcaction-diffusion-advection equations: 

^ + AfVu • Vu = F(u; /i) + DV^u, (2) 

where F is reaction kinetics, M represents a tensor of advection coefficients induced by endoplasmic fiow, and Z? is a 
diagonal matrix of diffusion constants. 

We assume that the system has a trivial homogeneous steady solution u = o, and it is stable when the bifucation 
parameter \s ^ < jic- Over the bifurcation point, /i > /ic, the trivial solution is unstable and local reaction kinetics 
becomes oscillatory. Reaction kinetics F(u; /i) can be expanded around the trivial solution as 

F(u; ^) = Lu + N2UU + iVsUuu -f • ■ • , 

^_ ^ _ f A N - ( \ at u-o 

Linearizing the system near the homogeneous steady state, we obtain equations for spatial Fourier components. 
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W{t; k) = {L- k^D)W{t; k), u{t, x) = / W{t; k)e'^''dk 



The stability of W = o is given by the eigenvalue problem for the eigenvalue A and eigenvector U, 

iL--k^D)\J ^ XIJ. 



(3) 



We obtain the eigenvalue A = (U*,(L - fc^Z?)U)/(U*, U) from Eq. (||). On the bifucation point, ReA = and 
d(ReX) / dk = are satisfied at some k — kc for the maximal eigenvalue(s). In the vicinity of the bifurcation point 
(/Lt, k) = (/ic, kc), the eigenvalue is 



/ ^ 9X ^ . I d'^X ^, , ^2 



(4) 



where the subscript c denotes the bifurcation point, and expansion coefficients are given as 



dX 
dk 



= ±ico 



UJc, Cg > 0. 



Hence the bifurcation problem has four generic cases: (i) kc = and lUc = 0; (ii) kc = and ojc 7^ (Hopf bifurcation) ; 
(iii) kc ^ and lOc = (Turing bifurcation); (iv) kc ^ and 7^ (travelling- wave type). According to the weakly 
nonlinear analysis of Eq. (||), the advection term is expected to modulate its bifurcation behaviour in the case of (iii) 
and (iv). 



III. ENVELOPE EQUATION 

Let us consider the envelope equation just after the bifurcation of travelling-wave type occurs, on a basis of weakly 
nonlinear analysis We denotes the bifurcation parameter by (/x — /ic) ~ e^, and suppose that ~ 0(e) 

(j = 1, . . . , A'^) in the vicinity of the bifurcation point. In the following analysis, the envelope equation is derived for 
one spatial dimension system of Eq. (^ with the single travelling wave. We introduce perturbation expansions and 
multiple scales, 

u eui + e^U2 + e''^U3 H , L - Lq + £^-^2 H , 

X = X — Cpt, ^ — e{x — Cgt), T = e'^t, (5) 

where Cp = LOc/kc is the phase velocity and Cg is the group velocity. Substitution of Eq. (^) into (^) yields perturbation 
equations for each order in e: 

0(e™) (io + i^—^ + Cp — )u„, =b,„, m= 1,2,3,... (6) 

where denotes the inhomogeneous term of the m-th order equation. 

For the first order equation in Eq. (^, the inhomogeneous term is bi = o. Then we have a solution 

ui = W^(^,T)e''^'<=^U-hc.c., 

where c.c. means complex conjugate. U is an eigenvector of an eigenvalue A = —iuc for the eigenvalue problem (|^) 
on the bifurcation point. 

For the second order equation, that is m = 2 in Eq. (|^), we expand the solution and inhomogeneous term by the 
phase (j) = kcX — Wct, 



Then the solvability conditions for U2 are 



(U*,b(+^))-(U*,b^^))=0, (7) 
5C 



ba = — —icg + 2ikcD)\J, hY ' = b^ 
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These conditions are obviously satisfied since 



.dx 



(u*,u) • 



Thus we must advance our calculation to the third order to obtain the envelope equation. 
For the third order equation in e, we expand U3 and ba as 



U3 = Eu(')e^'^ b3=5:b^"e''^ 



then the solvability conditions for U3, 



(U*,b(+^)) = (U*,b("^') = 0, (8) 



b(+i) = - + 2klWM{n'^i^\J + Uu^^)) - WL,\J 



-|M^piyA^3(UUU + UUU + UUU) - 2ikcD-^ 9^^^' 



b(-).b(-). 

Thus we obtain the complex Ginzburg-Landau (CGL) equation: 
dW o d'^W 

— ^c,W-C2\W\'W + c,^, c, = (U*,V,)/(U*,U), (9) 
Vi = L2U, 

+7V2{tJ, (io^'^V^ (iV^+^%U) } - TVa (uUU + UUU + UUu) , 

V3 = -(cg + 2ik,D){Ll^'''>)-\cg + 2ik,D)V, 

where ^ Lq - [Ik^fD + imujc, n!^^ ^ N2 + Ik^M, and {X, Y} ^ XY + YX. By means of the dispersion 



relation (|4|), ci = {dX/dfi)c and C3 = {l/2){d'^X/dk^)c are generically complex constants. 

The CGL equation describes the small amplitude dynamics of the system near the bifurcation point, that is, 
W slowly and slightly modulates travelling wave with the wave number kc and frequency loc- The coefficient of the 
nonlinear term, C2, depends on the advection term in Eq. (|^). Since the advection term has the form of the gradient 
of metabolic species, C2 is including the term of k^. It is known that C2 determines the amplitude of plane wave 
solutions and the nonlinear dispersion relation. The signature of Re(c2) decides the type of bifurcation at Re(ci) = 0: 
a supercritical bifurcation occurs for Re(c2) > 0, while subcritical one for Re(c2) < 0. In the latter case, we will need 
higher order terms as to the CGL equation. 

The CGL equation has various type of solutions with relation to the coefficients , but wc do not mention each of 
them. Because we have derived the CGL equation from Eq. (|^) which has no concrete form of reaction and advection 
terms, it is impossible to give the coefficients of the CGL equation explicitly. In the next section, we discuss the 
advection effect implied by the CGL equation ^ and compare our results with ones of reaction-diffusion models 
without the endoplasmic flow. 

We comment on the envelope equation of spatial two dimensions. Add a scaling of y-axis, 77 — ey, to Eq. and 
continue the similar calculation above, we obtain the CGL- type equation without rotational symmetry (see appendix). 



IV. CONCLUSIONS 



We have presented a reaction-diffusion-advection model of the Physarum plasmodium which makes rhythmic con- 
traction. The local contraction determines intracellular transport of endoplasmic sol, and this motivates the crawling 
behaviour of the plasmodium. By means of the weakly nonlinear analysis of the model, we obtain the following 
results: (1) an advection effect arises to the slow dynamics when the travelling-wave type bifurcation occurs; (2) 
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the envelope equation has the form of the CGL equation; (3) the advection terms affect the nonhnear term of the 
CGL equation through k^. In weakly nonlinear region of a reaction-diffusion- advection model, (2) implies that we 
can observe similar dynamical behaviour to conventional reaction-diffusion systems without flow. Thus it is expected 
that some experimental results of the rhythmic contraction in the plasmodium have been illustrated by such simple 
reaction-diffusion systems. The last result (3), however, shows possibility that the nonlinear effect near the bifurcation 
point stems from not only reaction kinetics but also the self-determined flow. In the remainder of this section, we 
attempt to discuss the effect of the self-determined flow on physiological behaviour. 

In general, advection can play an effective role for pattern formation, as the flow of matter often causes instabilities 
of hydrodynamical systems |2^,Q. For example, Rovinsky and Menzinger have shown that a differential flow of 
chemical species induces instabilities of the homogeneous steady state, and leads to travelling wave pattern nevertheless 
without diffusions. In this case, advection terms have a crucial effect on linear dispersion relation of the wave 
propagation. Contrary to this, the advection studied in the present paper, the self-determined flow, has no effect on 
the linear stability, but modulates nonlinear dispersion relation. 

In weakly nonlinear region, the self-determined advection causes the strong dependency of pattern formation on 
the critical wave number, kc. The affected nonlinear term of the CGL equation is known to be closely related with 
bifurcation and size of amplitude. Because phase difference of the metabolic oscillation between neighbours tends 
to become larger as kc is larger, kc is regarded as an indicator that an oscillator is in step with its neighbouring 
oscillators. This point of view implies that phase difference in local plays an important role for pattern formation 
such as amoeboid behaviours in the Physarum plasmodium. 

Real patterns observed by measuring thickness oscillation in the plasmodium show that the wave number is always 
about one or two in free locomotion under culture conditions. However, larger wave number more than ten appears 
transiently in relation to changes in behaviours induced by an appropriate stimulation. This phenomenon may relate 
to the /cc-dependency. Further analysis of the kc-eSect would be given if the model equations are concretely specified. 
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APPENDIX: DERIVATION OF THE ENVELOPE EQUATION IN TWO DIMENSIONS 

Before studying the envelope equation in two dimensions, we comment on counter propagating waves. For the 
bifurcation of travelling-wave type, counter propagating waves are possible in one spatial dimension, although we 



derived the envelope equation (CGL equation) for the single travelling wave in section III. The envelope dynamics of 



the counter propagating waves is the coupled equations of envelopes W+ and W- of a linearlized solution 

where {Lq — k'^D)XJ± — ±iWcU±. The time evolution of W± is governed by the dynamics like CGL equation, but is 
modulated by the interaction terms |VF:^pM^±. 

In the case of two spatial dimensions, multimode travelling waves satisfying |kj| = kc, 



(Wj-i-lJ+ exp [i(kj • X + LUct)] + Wj-IJ- exp [i (k^ • x — LUct)] 



are possible because modes corresponding to an anulus of wave vectors, |k| — kc, are neutrally stable on the bifurcation 
point. Although we can obtain the multimode envelope equations, it is not clear that the dynamics of them has 
structural stability. Such a problem has a relation to pattern selection, and we need more precise analysis on the 
bifurcation with symmetry. Hereafter we discuss the envelope equation for a single travelling wave of Eq. (|2|) in two 
spatial dimensions. We assume that a set of mode near the single plane wave travelling along x-axis, exp[i{kcX — uJct)], 
is dominant across the bifurcation point. 

Introducing perturbation expansions and multiple scales, 

u ^ eui + e^U2 + e^U3 H , L ^ Lq + €^L2 H , 

X = X — Cpt, ^ = e{x — Cgt), ?7 = ey, t ~ e^t, (Al) 
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and substituting Eq. (Al) into Eq. (g), we obtain perturbation equations, 

0(e") (io + ^^+Cp^)u„ =b„, m= 1,2,3,... (A2) 

The inhomogeneous term of the first order equation is vanishing, bi = o. For the second and third order equations, 
the inhomogeneous terms are 

dui dui dui d^ui 

"2 — —Cg—^ + M ^^-^TT^ — iV2UiUi — 2D 



dX dX dXdC 

dui du2 /9u2c'ui 9ui 9u2 9ui 9ui 9ui(9ui\ 

( <9^U2 (9^ui 9^ui 
-L2U1 - iV2(u2Ui + U1U2) - A^sUiUiUi - Z?! 2 



We write a sohition of the first order equations with the slowly varying envelope W{S,, rj, t), 

ui = We:''''^V + c.c. 

Since the solvability conditions for the second order equation are satisfied by the definition of Cg, we advance our 
calculation to the third order: 

(U*,b(+^)) = (U*,b(-^))=0, (A3) 
b(+i) = - + 2fc,^iyAf(uf U + Uu(^)) - WL,V 

-T4^iV2(u^°^U + Uu^°^) - M^iV2(u^^^U + IJu^^^) 

i^ii'"^^ r)'^W\ 
-\W\^Wm{VlJiJ + UUU + UUU) - 2tk,D^ - + ^)DU, 



From the solvability conditions (A3), we get the envelope dynamics, 

— = ciW^ - c,\W\'W + (C3 + C4)^ + C4^, (A4) 

where coefficients ci, C2, C3 are the same in Eq. ^ and C4 — —icg/2kc. The coefficient C4 implies dipersion waves 
which are observed in the nonlinear Schrodinger equation. 
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